home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_2 / a2_10.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  4.5 KB  |  189 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 2.10 (Newton-Raphson Method in 2-Dimensions).
  9. % Section    2.7,    Newton's Method for Systems, Page 116
  10. echo on; clc; format long; hold off; clear
  11. % This program implements the Newton-Raphson method
  12.  
  13. % for solving    0 = f1(x,y)  and  0 = f2(x,y).
  14.  
  15. %    Use the vector notation X = (x y).
  16.  
  17. %    Define and store the functions f1(X) and f2(X) and
  18.  
  19. %    define and store the functions  F(X) and  J(X) in the
  20.  
  21. %    M-files  f1.m   f2.m   F.m   and   J.m   respectively.
  22.  
  23. pause    % Press any key to continue.
  24.  
  25. clc;
  26. % function z = f1(x,y)
  27. % z = x.^2 - 2.*x - y + 0.5;
  28.  
  29. % function z = f2(x,y)
  30. % z = x.^2 + 4.*y.^2 - 4;
  31.  
  32. % function Z = F(X)
  33. % x = X(1); y = X(2);
  34. % Z = [f1(x,y);  f2(x,y)];
  35.  
  36. % function W = J(X)
  37. % x = X(1); y = X(2);
  38. % W = [(2.*x - 2)  (-1);
  39. %      (2.*x)      (8.*y)];
  40.  
  41. delete f1.m
  42. diary f1.m; disp('function z = f1(x,y)');...
  43.             disp('z = x.^2 - 2.*x - y + 0.5;');...
  44. diary off;
  45.  
  46. delete f2.m
  47. diary f2.m; disp('function z = f2(x,y)');...
  48.             disp('z = x.^2 + 4.*y.^2 - 4;');...
  49. diary off;
  50.  
  51. delete F.m
  52. diary F.m; disp('function Z = F(X)');...
  53.            disp('x = X(1); y = X(2);');...
  54.            disp('Z = [f1(x,y); f2(x,y)];');...
  55. diary off;
  56.  
  57. delete J.m
  58. diary J.m; disp('function W = J(X)');...
  59.            disp('x = X(1); y = X(2);');...
  60.            disp('W = [(2.*x - 2)  (-1);');...
  61.            disp('     (2.*x)      (8.*y)];');...
  62. diary off;
  63.  
  64. % Remark. f1.m, f2.m, F.m, J.m and new2dim.m are used for Algorithm 2.10
  65. f1(0,0); f2(0,0); F([0 0]); J([0 0]); % Test for files f1.m f2.m F.m J.m
  66. pause      % Press any key to graph 0 = f1(x,y) and 0 = f2(x,y).
  67.  
  68. clc; clg;
  69. a = -2.0;
  70. b = 2.5;
  71. c = -1.0;
  72. d = 1.5;
  73. hx = (b-a)/31;
  74. hy = (d-c)/31;
  75. [X Y] = meshdom(a:hx:b, c:hy:d);...
  76. Z = f1(X,Y);...
  77. W = f2(X,Y);...
  78. contour(Z,[0 0],a:hx:b,c:hy:d,'-g');...
  79. hold on;...
  80. contour(W,[0 0],a:hx:b,c:hy:d,'-g');...
  81. title('Implicit plot of 0 = f1(x,y) and 0 = f2(x,y).');...
  82. grid;...
  83. hold off;...
  84. shg; pause    % Press any key to perform Newton-Raphson iteration.
  85.  
  86. clc;
  87. %    Place the starting value in  [p0 q0]
  88.  
  89. % Place the tolerance for P in delta
  90.  
  91. % Place the tolerance for F(P) in epsilon
  92.  
  93. %    Place the maximum number of iterations in  max1
  94.  
  95. p0 = 2.0;
  96. q0 = 0.25;
  97. P0 = [p0 q0];
  98. delta = 1e-12;
  99. epsilon = 1e-12;
  100. max1  = 40;
  101.  
  102. [P0,F0,err,P] = new2dim('F','J',P0,delta,epsilon,max1);
  103.  
  104. pause    % Press any key for the list of iterations.
  105.  
  106. clg; clc;
  107. a = 1.86;
  108. b = 2.02;
  109. c = 0.24;
  110. d = 0.34;
  111. hx = (b-a)/20;
  112. hy = (d-c)/20;
  113. [X Y] = meshdom(a:hx:b, c:hy:d);...
  114. Z = f1(X,Y);...
  115. W = f2(X,Y);...
  116. contour(Z, [0 0],a:hx:b,c:hy:d,'-g');...
  117. grid;...
  118. hold on;...
  119. contour(W, [0 0],a:hx:b,c:hy:d,'-g');...
  120. X0 = P(:,1);...
  121. Y0 = P(:,2);...
  122. plot(X0,Y0,'or');...
  123. title('Graphical presentation of the Newton-Raphson iteration.');...
  124. hold off;...
  125. shg; pause    % Press any key to continue.
  126.  
  127. Mx1 = 'Iterations for Newton`s method.';
  128. Mx2 = '     p(k)               q(k)';
  129. Mx3 = 'The solution is:';
  130. Mx4 = 'The error estimate for P is ';
  131. Mx5 = num2str(err);
  132. Mx6 = [Mx4,'[± ',Mx5,', ± ',Mx5,']'];
  133. clc,echo off, diary output,...
  134. disp(''), disp(Mx1),disp(''), disp(Mx2),disp(P),...
  135. disp('Iteration converged to the root.'),...
  136. disp(''),disp(Mx3),disp(''),disp('P = '),...
  137. disp(P0),disp('F(P) = '),disp(F0),...
  138. disp(''),disp([Mx6]),diary off,echo on
  139.  
  140. pause    % Press any key to perform Newton-Raphson iteration.
  141.  
  142. clc;
  143. %    Place the starting value in  [p0 q0]
  144.  
  145. % Place the tolerance for P in delta
  146.  
  147. % Place the tolerance for F(P) in epsilon
  148.  
  149. %    Place the maximum number of iterations in  max1
  150.  
  151. p0 = -0.1;
  152. q0 =  0.9;
  153. P0 = [p0 q0];
  154. delta = 1e-12;
  155. epsilon = 1e-12;
  156. max1  = 40;
  157.  
  158. [P0,F0,err,P] = new2dim('F','J',P0,delta,epsilon,max1);
  159.  
  160. pause    % Press any key for the list of iterations.
  161.  
  162. clg; clc;
  163. a = -0.26;
  164. b = -0.08;
  165. c =  0.88;
  166. d =  1.02;
  167. hx = (b-a)/20;
  168. hy = (d-c)/20;
  169. [X Y] = meshdom(a:hx:b, c:hy:d);...
  170. Z = f1(X,Y);...
  171. W = f2(X,Y);...
  172. contour(Z,[0 0],a:hx:b,c:hy:d,'-g');...
  173. grid;...
  174. hold on;...
  175. contour(W,[0 0],a:hx:b,c:hy:d,'-g');...
  176. X0 = P(:,1);...
  177. Y0 = P(:,2);...
  178. plot(X0,Y0,'or');...
  179. title('Graphical presentation of the Newton-Raphson iteration.');...
  180. hold off;...
  181. shg; pause    % Press any key to continue.
  182.  
  183. clc,echo off, diary output,...
  184. disp(''), disp(Mx1),disp(''), disp(Mx2),disp(P),...
  185. disp('Iteration converged to the root.'),...
  186. disp(''),disp(Mx3),disp(''),disp('P = '),...
  187. disp(P0),disp('F(P) = '),disp(F0),...
  188. disp(''),disp([Mx6]),diary off,echo on
  189.